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ABSTRACT 

We present a survey of 41 Kepler Objects of Interest (KOIs) for exomoons using Bayesian photody¬ 
namics, more than tripling the number of KOIs surveyed with this technique. We find no compelling 
evidence for exomoons although thirteen KOIs yield spurious detections driven by instrumental arti¬ 
facts, stellar activity and/or perturbations from unseen bodies. Regarding the latter, we find seven 
KOIs exhibiting > 5 cr evidence of transit timing variations, including the “mega-Earth” Kepler-lOc, 
likely indicating an additional planet in that system. We exploit the moderately large sample of 
57 unique KOIs surveyed to date to infer several useful statistics. For example, although there is a 
diverse range in sensitivities, we find that we are sensitive to Pluto-Charon mass-ratio systems for 
~ 40% of KOIs studied and Earth-Moon mass-ratios for 1 in 8 cases. In terms of absolute mass, our 
limits probe down to 1.7 Ganymede masses, with a sensitivity to Earth-mass moons for 1 in 3 cases 
studied and to the smallest moons capable of sustaining an Earth-like atmosphere (0.3 M 0 ) for 1 in 4. 
Despite the lack of positive detections to date, we caution against drawing conclusions yet, since our 
most interesting objects remain under analysis. Finally, we point out that had we searched for the 
photometric transit signals of exomoons alone, rather than using photodynamics, we estimate that 1 
in 4 KOIs would have erroneously been concluded to harbor exomoons due to residual time correlated 
noise in the Kepler data, posing a serious problem for alternative methods. 

Subject headings: techniques: photometric — planetary systems 


1. INTRODUCTION 

The occurrence rate, and even the existence, of moons 
orbiting extrasolar planets is presently unknown. Their 
discovery would provide, for the first time, information 
on the physical processes governing satellite formation 
and evolution outside the Solar System. Additionally, 
even null detections can guide theory in this respect, 
when carefully interpreted in terms of upper limits using 
the appropriate statistical techniques. Since it is a-priori 
unknown whether exomoons are sufficiently large to be 
detected using current instrumentation, this realization 
motivates a survey program to be designed such that 
statistically rigorous upper limits are derived for null de¬ 
tections. Whilst this strategy necessitates substantially 
greater computational resources, a survey adopting this 
philosophy has a guaranteed science product even if no 
confirmed objects are ever found, which in this case is 
the occurrence of exomoons as a function of the planet 
and satellite properties. 

The “Hunt for Exomoons with Kepler” (HEK) project 
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(Kipping et al. 2012) is a survey for exomoons around 
transiting exoplanets which follows this ethos. Our sur¬ 
vey uses Bayesian inference to derive rigorous upper lim¬ 
its in the cases of null detections and employs a suite of 
vetting strategies to interrogate candidate signals, honed 
by our past experience of conducting this survey. In pre¬ 
vious works, we surveyed 18 Kepler Objects of Interest 
(KOIs) for evidence of exomoons including a sample of 
M-dwarfhost stars (Kipping ct al. 2014a), habitable-zone 
KOIs (Kipping et al. 2013b), dynamically active systems 
(Ncsvorny ct al. 2012) and single KOIs (Kipping ct al. 
2013a). In each case, we searched for both the transits of 
the exomoon and the perturbations induced by its mass 
with a self-consistent, analytic, dynamical model known 
as LUNA (Kipping 2011). Correlated noise structure is 
a particularly frequent source of exomoon transit false- 
positives, as outlined in the cautionary tale of Kepler- 
90g.01 (Kipping et al. 2015). 

In this paper, we dramatically increase the number of 
planetary candidates surveyed by HEK, presenting an 
analysis of 41 KOIs. Since two of these were already 
studied by previous surveys (but with fewer data), this 
paper brings the total number of unique KOIs surveyed 
for exomoons by HEK to 57. The paper is split as fol¬ 
lows. In §2, we describe our target selection strategy 
for these 41 KOIs via three distinct pathways. In §3, 
the data handling, modeling and inference methods are 
described, including the vetting and classification strate¬ 
gies we employ. In §4, we present the results of our sur¬ 
vey, including refined planet parameters and an analysis 
of transit timing variations. We conclude in §5 with a 
discussion focussed on some statistical results from this 
sample and the HEK surveys to date, including empirical 
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sensitivity limits and computational demands. 

2. TARGET SELECTION 

2.1. Overview 

As outlined in Kipping ct al. (2012), the HEK project 
ntilizes three complementary methods of selecting tar¬ 
gets for analysis: 

■ TSA: Target Selection Automatic 

■ TSV: Target Selection Visual 

■ TSO: Target Selection via Opportunities 

In previous papers, we limited our samples to targets 
selected using just one of these methods. In this work, 
thanks to the larger number of KOIs under consideration, 
we selected mnltiple targets from each category. Specif¬ 
ically, our sample consists of 16 KOIs selected via TSA, 
16 via TSV and 9 via TSO, as described below. The full 
list of KOIs considered here is listed in Table 1. 

2.2. Target Selection Automatic (TSA) 

The TSA algorithm is unchanged from that described 
in Kipping et al. (2014a) and references therein. In sum¬ 
mary, the algorithm assigns an estimated mass and tidal 
properties (specifically k 2 p/Q) to each KOI based on 
the size from empirical and theoretical relations. As¬ 
suming loss through tidal evolution, we are then able to 
compute the maximum moon mass that can survive for 
5Gyr using the expressions of Barnes & O’Brien (2002). 
This maximum moon mass is then converted into a ra¬ 
dius using an appropriate mass-radius relation and finally 
into a signal-to-noise ratio (SNR) using the Kepler Com¬ 
bined Differential Photometric Precision (CDPP) statis¬ 
tics (Christiansen ct al. 2012). 

Target KOIs were selected from the NASA Exoplanet 
Archive list of KOIs, excluding objects flagged as false- 
positives. In this work, we elect to restrict the TSA tar¬ 
gets to those KOIs with radii Rp > 2R^, insolation 
SeS < 6.925 50 (inner edge of the habitable-zone defined 
by Zsom et al. 2013) and SNR > 4/7075 = 4.62. The 
70.75 scaling is based on the fact TSA computes SNR 
assuming one uses the entire time series, but in reality we 
only use 75% of the data in our fits for reasons described 
in §3. 

In total, 25 targets satisfied these criteria, of which 
three have already been studied in previous surveys 
(Kepler-22b, Kipping ct al. 2013b; KOI-1596.02, Kip¬ 
ping ct al. 2014a; KOI-1876.01, Kipping ct al. 2013a). 
Since significantly more data is now available for KOI- 
1876.01 than available for the analysis in Kipping et 
al. (2013a), this target was kept in our TSA sample. 
KOI-1871.01 and KOI-2762.01 were removed from this 
subset because they were flagged as candidates by the 
TSV process described later in §2.3. Additionally, KOI- 
1174.01, KOI-490.02 and KOI-1274.01 (also known as 
Kepler-421b; Kipping ct al. 2014b) have fewer than three 
visible transits and were rejected. Finally, KOI-902.01 
was reported by klazch ct al. (2013) to display high- 
amplitude (>100 mins) and long-period (> 3 years) tran¬ 
sit timing variations (TTVs), which are uncharacteristic 
of exomoons (Kipping 2009a), and so was also rejected. 
This resulted in the 16 KOIs listed in Table 1. 


2.3. Target Selection Visual (TSV) 

In previous HEK papers, the only target picked for 
analysis using our TSV process was KOI-0872.01, also 
known as Kepler-46b (see the Supporting Online Mate¬ 
rial of Nesvorn/" et al. 2012). This target displayed in¬ 
transit flux increases characteristic of star-planet-moon 
mutual events (as well as large TTVs), but these were ul¬ 
timately attributed to star-spot crossings (e.g. see Bcky 
ct al. 2014) rather than due to an exomoon. 

In the earlier work, TSV was conducted using sev¬ 
eral publicly available light curve viewers. However, we 
quickly found that these viewers did not meet our needs 
for performing large scale visual surveys. In a typical 
survey, four to five thousand KOI transit signals need to 
be examined for exomoon signals across several hundred 
distinct planetary candidates. To effectively handle such 
a large volume of data in a reasonable time frame requires 
a system with advanced navigation and signal processing 
features, which the public tools did not provide. 

To address these needs, we developed a custom de¬ 
signed software package known as LcTools. This package 
consists of two main components: LeCenerator for gen¬ 
erating normalized light curve files in bulk, and LeViewer 
for examining light curve files and recording signals. This 
Windows-based software is freely available and can be 
obtained from the author (A. Schmitt) upon request. We 
direct those interested to the online description for tech¬ 
nical details^^, but we briefly point out some relevant 
features here: 

■ Bin short-cadence (SC) data to a user specified ca¬ 
dence. 

■ Independently pan and scale the time and flux axes. 

■ Locate each instance of a periodic signal and automat¬ 
ically scale the signal for optimized viewing. 

■ Automatically highlight known transit signals accord¬ 
ing to the NASA Exoplanet Archive ephemerides (e.g. 
see Figure 1). 

■ Automatically highlight mutual planetary transit sig¬ 
nals (the simultaneous transit of two or more planets). 

■ Record and highlight user-defined signals (e.g. due to 
exomoons), as shown in Figure 1. 

In total, 1289 unique KOIs were examined for exomoon 
signals across numerous surveys. Many KOIs were in¬ 
spected multiple times. KOIs included both TSAs and 
non-TSAs. The latter group was largely composed of 
promising targets meeting the selection criteria below: 

► Semi-major axis; a > 0.1 AU (for reasons discussed in 
Namouni 2010). 

► Planet size; Rp > 3i?0. 

► Expected SNR for an exomoon satisfies SNR > 1. 

In all cases, the minimum transit duration required was 
four hours. From experience, we find that at least eight 

See http://www.exomoon.org for details 
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TABLE 1 

List of KOIs surveyed for exomoons in this work, and their basic parameters taken from the NASA Exoplanet Archive. 


KOI 

Pp[d] 

[-R©] 

SNR 

Multip. 

Teq [K] 

Teff [K] 

Kp 

Planet Name 

Validation Ref. 

16 via TSA 

KOI-0438.02 

52.7 

2.06 

6.77 

2 

298 

3985 

14.3 

Kepler-155c 

Rowe et al. (2014) 

KOI-0518.03 

247.4 

2.14 

9.74 

3 

212 

4590 

14.3 

Kepler-174d 

Rowe et al. (2014) 

KOI-0854.01 

56.1 

2.05 

5.73 

1 

252 

3593 

15.8 

- 

- 

KOI-0868.01 

236.0 

8.73 

7.13 

1 

172 

3822 

11.6 

- 

- 

KOI-1361.01 

59.9 

2.53 

7.38 

1 

298 

4014 

15.0 

Kepler-6 lb 

Ballard et al. (2013) 

KOI-1431.01 

345.2 

8.00 

4.65 

1 

272 

5628 

13.5 

- 

- 

KOI-1783.02 

284.0 

4.80 

4.79 

2 

321 

6215 

13.9 

- 

- 

KOI-1830.02 

198.7 

2.88 

6.03 

2 

257 

5080 

14.4 

- 

- 

KOI-1876.01 

82.5 

2.56 

7.78 

1 

276 

4237 

15.2 

- 

- 

KOI-2020.01 

111.0 

2.32 

10.17 

1 

260 

4441 

15.5 

- 

- 

KOI-2686.01 

211.0 

3.68 

7.21 

1 

233 

4628 

13.9 

- 

- 

KOI-2691.01 

97.5 

3.30 

6.35 

1 

313 

4728 

15.0 

- 

- 

KOI-3263.01 

76.9 

7.00 

5.03 

1 

211 

3587 

15.9 

- 

- 

KOI-4005.01 

178.1 

2.22 

5.55 

1 

316 

5592 

14.6 

Kepler-439b 

Torres et al. (2015) 

KOI-4036.01 

168.8 

2.24 

5.93 

1 

279 

4893 

14.1 

- 

- 

KOI-4054.01 

169.1 

2.21 

5.58 

1 

308 

5380 

14.6 

- 

- 

16 via TSV 

KOI-0092.01 

65.7 

3.07 

8.70 

1 

503 

5952 

11.7 

- 

- 

KOI-0112.01 

51.1 

3.04 

6.74 

2 

552 

5803 

12.8 

- 

- 

KOI-0209.01 

50.8 

9.16 

5.29 

2 

636 

6466 

14.3 

Kepler-117b 

Rowe et al. (2014) 

KOI-0276.01 

41.7 

2.71 

5.86 

1 

629 

5982 

11.9 

- 

- 

KOI-0308.01 

35.6 

2.88 

2.28 

1 

672 

6210 

12.4 

- 

- 

KOI-0374.01 

172.7 

3.14 

9.26 

1 

389 

5839 

12.2 

- 

- 

KOI-0398.01 

51.8 

8.77 

2.98 

3 

459 

5271 

15.3 

- 

- 

KOI-0458.01 

53.7 

8.16 

2.65 

1 

558 

5833 

14.7 

- 

- 

KOI-0847.01 

80.9 

4.83 

1.13 

1 

413 

5665 

15.2 

- 

- 

KOI-1535.01 

70.7 

2.13 

7.67 

1 

539 

6174 

13.0 

- 

- 

KOI-1726.01 

45.0 

2.01 

9.63 

1 

405 

4684 

13.1 

- 

- 

KOI-1808.01 

89.2 

3.82 

5.36 

1 

504 

6277 

12.5 

- 

- 

KOI-1871.01 

92.7 

2.25 

7.34 

2 

295 

4528 

14.9 

- 

- 

KOI-2065.01 

80.2 

3.60 

4.33 

1 

437 

5685 

14.3 

- 

- 

KOI-2762.01 

133.0 

2.15 

4.69 

1 

266 

4530 

15.0 

- 

- 

KOI-5284.01 

389.3 

6.79 

2.10 

1 

267 

5731 

14.6 

- 

- 

9 via TSO 

KOI-0072.02 

45.3 

2.30 

12.44 

2 

554 

5627 

11.0 

Kepler-lOc 

Dumusque et al. (2014) 

KOI-0245.01 

39.8 

1.88 

49.27 

3 

456 

5417 

9.7 

Kepler-37d 

Barclay et al. (2013) 

KOI-0386.02 

76.7 

2.89 

4.15 

2 

512 

6197 

13.8 

Kepler-146c 

Rowe et al. (2014) 

KOI-0518.02 

44.0 

1.43 

10.37 

3 

377 

4590 

14.3 

Kepler-174c 

Rowe et al. (2014) 

KOI-0722.01 

46.4 

2.85 

3.12 

1 

635 

6343 

13.5 

- 

- 

KOI-1783.01 

134.5 

7.93 

6.06 

2 

412 

6215 

13.9 

- 

- 

KOI-2358.01 

56.5 

2.04 

6.39 

1 

583 

6412 

13.5 

- 

- 

KOI-3663.01 

282.5 

11.28 

5.20 

1 

328 

6007 

12.6 

PH-2b 

Wang et al. (2013) 

KOI-3681.01 

217.8 

22.02 

2.58 

2 

497 

6382 

11.7 

- 

- 


long-cadence (LC) data points are necessary to reliably 
detect exomoon signals visually. 

For a given host star, all available Presearch Data Con¬ 
ditioning (PDCSAP) data (Stumpe ct al. 2012, Smith et 
al. 2012) were obtained from the Barbara A. Mikulski 
Archive for Space Telescopes^^ (MAST). The data were 
normalized and merged into a single light curve file. Both 
the LC and SC data were used, with SC given priority. 
To mitigate the flux scatter and improve signal detection, 
SC data were binned to a lower data rate - typically 6, 
10, or 12 data points per hour. 

For each KOI transit signal, a visual search was made 
for moon-like signals within the marked KOI signal re¬ 
gion (see Figure 1 for an example). The KOI signal region 
is defined as the extended Hill timescale per the equation 
given in §2.2 of Kipping ct al. (2013b). Upon finding a 

https://archive.stsci.edu/kepler 


moon-like signal conforming to the general shape char¬ 
acteristics predicted by LUNA (Kipping 2011), the signal 
was manually recorded in LcViewer and a screen snap¬ 
shot taken. 

A total of 2928 moon-like signals were recorded; 923 
from the TSA group and 2005 from the non-TSA group. 
From this set, the highest quality signals were visually 
identified. For the current paper, 16 KOIs were found 
having at least three high quality exomoon signals (see 
Table 1). Note that this represents just a fraction of the 
KOIs with high-quality signals found. 

2.4. Target Selection via Opportunities (TSO) 

In addition to the 32 KOIs selected using the TSA and 
TSV methods, we picked 9 KOIs as “targets of opportu¬ 
nity” via TSO. This process is loosely defined as iden¬ 
tifying systems for which we have reason to suspect an 
elevated probability of a successful detection or that the 
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Fig. 1.— Main application window for LcViewer. Two signals are shown in this example: a planetary transit signal of KOI-0308.01 
highlighted in orange and a user defined moon-like signal highlighted in green. The white vertical lines on either side of the light curve 
indicate the region of visual inspection, which is equal to the estimated Hill timescale. 


system has some unusual properties making an analysis 
of greater scientific value. As an example from previous 
surveys, in Kipping ct al. (2013b) we selected Kepler-22b 
using TSO, since it was claimed to be the first habitable- 
zone transiting planet discovered (Borucki et al. 2012). 
Below, we briefly describe the reasons why each KOI was 
selected using TSO. 

KOl-0722.01 was previously surveyed in Kipping et al. 
(2013a) using quarters Ql-9 of the Kepler data, but is re¬ 
analyzed here using Ql-17 due to the unusual fits found 
previously. 

KOl-0072.01, also known as Kepler-lOc, is a confirmed 
transiting planet with a mass measured using radial ve¬ 
locities (Dumusquc ct al. 2014). Argued by the authors 
to be a solid Neptune-mass planet, Kepler-lOc may repre¬ 
sent a scaled up version of the Earth leading us to wonder 
if there is a scaled up version of the Moon in orbit too. 

KOl-3663.01, also known as PH-2b, is a validated tran¬ 
siting planet (Wang ct al. 2013) found by the Planet 
Hunters project (Fischer ct al. 2012). The discoverers 
argue that the planet is a Jupiter-sized planet in the 
habitable-zone of its host star, making it an ideal target 
for an exomoon search. 

KOl-1783.01 and KOI-0518.02 (also known as Kepler- 


174c; Rowe et al. 2014) are additional KOls to two 
targets already selected in the TSA part of this work, 
namely KOl-1783.02 and KOI-0518.03 (also known as 
Kepler-174d; Rowe ct al. 2014). Moreover, both of these 
have relatively high SNR expected due to putative exo¬ 
moons and so were selected via TSO. 

KOI-0245.01, also known as Kepler-37d, is a validated 
transiting planet and the largest and widest orbit planet 
in a three-planet system (Barclay ct al. 2013). The pho¬ 
tometric precision achieved for Kepler-37 is exceptional, 
thanks to the star’s brightness (9.7 in Kepler's bandpass) 
and low intrinsic noise, allowing for the discovery of the 
Moon-sized planet Kepler-37b. For this reason, we added 
this target through TSO. 

KOI-3681.01 was surveyed for exomoons in an inde¬ 
pendent study (Lewis 2014), but only using the transit 
method rather than the full photodynamic modeling em¬ 
ployed by HEK. Nevertheless, we added this target to 
exploit the unique opportunity to compare our results. 

Finally, KOl-0386.02 and KOl-2358.01 were selected 
as somewhat “wildcard” targets due to their moderate 
orbital periods providing plenty of transits (76.7 d and 
56.5 d respectively) and yet TSA still predicting reason¬ 
ably good SNR values of 4.15 and 6.39 respectively. 
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3. ANALYSIS METHODOLOGY 

3.1. Data Handling 

Photometric time series for the 41 KOIs studied in this 
work were obtained from the MAST archive. We down¬ 
loaded all available quarters from QO-17 for each KOI 
and preferentially use the SC over LC data, where avail¬ 
able. Only data with SAP QUALITY values equal to 
zero were utilized in this work, thereby rejecting data 
likely corrupted by known effects. 

As with previous papers, we only use the Simple Aper¬ 
ture Photometry (SAP) processed data in our analysis, 
rather than the Presearch Data Conditioning (PDCSAP) 
time series (Stumpe ct al. 2012, Smith et al. 2012). The 
PDCSAP data is optimized to remove common mode in¬ 
strumental effects in order to retain stellar astrophysics 
such as rotational modulations. However, in our work, 
retaining the stellar astrophysics is not required and crit¬ 
ically the act of cleaning common modes will remove fre¬ 
quencies across the Fourier domain that may potentially 
distort, remove or even spuriously introduce the very low 
amplitude exomoon-like signals we seek here. 

Transit light curves are detrended on an epoch-by¬ 
epoch basis. Transits are initially located using the 
NASA Exoplanet Archive reported ephemeris and we 
consider a window of ±0.5P around each event. Dis¬ 
continuities in the time series, due to pointing tweaks 
for example, are visually identified and in such cases our 
epoch window is shortened up to the nearest such event 
surrounding each transit i.e. we make no attempt to 
correct such effects. Long-term trends, due to focus drift 
for example, are removed using the CoFiAM algorithm de¬ 
scribed in Kipping et al. (2013a, 2014a). This algorithm 
is essentially a low-cut, high-pass filter removing all long¬ 
term trends with a threshold frequency set to three times 
the NASA Exoplanet Archive reported transit duration. 
As described in our previous papers, CoFiAM ensures the 
transit light curve of the planet and any moons are not 
distorted, whilst removing long-term trends and mini¬ 
mizing the autocorrelation on a 30 minute timescale. It is 
important to stress that CoF iAM does not remove high fre¬ 
quency noise, but our philosophy here is that whilst one 
can consider alternative methods to accomplish this feat, 
it is unclear whether such methods can retain the guar¬ 
antee of CoFiAM that the transit shape is undisturbed. 

3.2. Forward Modeling 

As with previous HEK papers, the forward model that 
is compared to the observations is the LUNA algorithm, 
described in Kipping (2011). At every time stamp, LUNA 
calculates the sky-projected positions of the planet and 
moon relative to the star, using a nested Keplerian body 
appropriate for a hierarchical three-body problem (see 
Appendix of Kipping 2010b for proof). Based on these 
positions, the amount of star light blocked by the planet 
and the moon are calculated analytically using the Fewcll 
(2006) and Kipping (2011) solutions. Limb darkening 
of the star is appropriately accounted for, provided the 
radius of the satellite is small relative to the star i.e. 
Rs R* (this is not required for the planet). In 
this work, we use the quadratic limb darkening formu¬ 
lation described in Kipping (2013) and account for long- 
cadence smearing using the re-sampling treatment de¬ 
scribed in Kipping (2010a) (A^resam = 30). Quarter-to- 


quarter contaminating light factors are taken from the 
fits file (“CROWDSAP” header keyword) and accounted 
for following the method described in Kipping & Tinetti 
( 2010 ). 

Recently, there have been several proposals for find¬ 
ing exomoons using alternative methods to photodynam¬ 
ics, such as excess noise around transits (“scatter-peak”; 
Simon ct al. 2012), high frequency TTVs (Szabo ct al. 

2013) and phase-stacking transits (“OSE”; Heller ct al. 

2014) , and we briefly comment on why we still favor the 
photodynamical technique. 

The alternative methods aim to be computationally 
cheaper than photodynamical modeling, at the expense 
of approximating the complex motions involved. How¬ 
ever, it is important to note that photodynamical model¬ 
ing accounts for all known photometric effects caused by 
a transiting exomoon, e.g., TTV (Sartorctti & Schnei¬ 
der 1999), TDV-V (Kipping 2009a), TDV-TIP (Kip¬ 
ping 2009b), plus even previously unknown effects such 
ingress/egress asymmetries (Kipping 2011). Therefore, 
these simpler methods can all be understood to be ap¬ 
proximations to the full photodynamical description and 
by definition model less of the information content im¬ 
parted into the light curve by the presence of exomoon(s). 
With less information, they are all guaranteed to not 
only be less sensitive in terms of signal-to-noise, but also 
critically lacking the powerful vetting techniques avail¬ 
able with photodynamics, such as dynamical weighing of 
the parent planet (Kipping 2010b). We also stress that 
our implementation of photodynamics models the entire 
Kepler time series, i.e. every transit is modeled simul¬ 
taneously and self-consistently. Naturally then, there is 
no more signal-to-noise to be gained by phase-folding or 
binning techniques (e.g. Simon et al. 2012; Heller et al. 
2014), which will actually only serve to further degrade 
the information content. For these reasons, we argue 
that photodynamics remains the most robust method of 
searching for exomoons. 

However, the alternative methods are typically much 
faster and may be useful for identifying interesting sys¬ 
tems worthy of more detailed analysis with photodynam¬ 
ics i.e. as a prioritization tool. The HEK project cur¬ 
rently adopts this position in its application of the Target 
Selection Visual process discussed in §2.3. 

3.3. Bayesian Inferenee & Model Selection 

For reasons discussed in §1, the HEK project places 
great importance on the use of rigorous statistical infer¬ 
ence techniques in our survey. There are several unique 
needs of an exomoon survey that shape our choice of 
methodology: 

■ The parameter space is complex, highly multi¬ 
modal and features curved degeneracies, posing a 
major challenge to many conventional regression 
techniques. 

■ In order to derive an upper limit on the satellite-to- 
planet mass ratio, (Ms/Mp), one must marginalize 
over the other parameters, such as the moon’s semi¬ 
major axis and the three-dimensional architecture 
of the orbits. 

■ A planet-Fmoon model has significantly greater 
freedom than a planet-only case guaranteeing a 
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higher maximum likelihood, thereby necessitating 
the use of Bayesian model selection. 

It should be recognized that satisfying these needs can¬ 
not be accomplished without high computational cost. 
One of the most efficient algorithms that satisfies these 
requirements is MultiNest (Fcroz ct al. 2008, 2009), 
based on nested sampling (Skilling 2004), and this is the 
method employed in this work and in our previous HEK 
papers. Using MultiNest, multiple objects (in our case 
multiple exomoons) are identifiable as separate modes 
(Fcroz et al. 2008, 2013), unless the moons strongly inter¬ 
act. This important feature, utilized since the inaugural 
HEK paper (Kipping ct al. 2012), refutes the claim made 
in Heller et al. (2014) that OSE is the first method per¬ 
mitting the characterization of multiple moon systems. 

For each KOI, we regress at least two models, or hy¬ 
potheses, to the data. The null hypothesis is "Hp, a 
planet-only case (forward model is that of Mandcl & Agol 
2002), which is compared to the planet-l-satellite hypoth¬ 
esis, T-Ls (forward model is LUNA). The parameterization 
of these models is the same as that described in Kipping 
et al. (2014a) and the priors adopted are listed in Ta¬ 
ble 2. For both models, a Gaussian likelihood function 
is adopted in our fits and exomoons are assumed to have 
near-circular orbits. 

As in our previous paper Kipping ct al. (2014a), we 
also have latent priors on the derived mean planetary and 
satellite densities, pp and ps (see Kipping 2010b for how 
these may be calculated). Specifically, we impose a max¬ 
imum density for both of 27.9 gcm“^ (see Kipping ct al. 
2012 for justification) and a minimum density for pp of 
0.08 gcm“^ (no such minimum is imposed on ps to allow 
null detections of Ms/Mp). The purpose of these latent 
priors is to restrict MultiNest to only the physically 
plausible range of solutions (thereby imposing detection 
criteria B3 automatically; see §3.4). One drawback of 
this decision is that a zero-radius moon (i.e. a null de¬ 
tection of a moon signal) is forbidden, since ps ^ oo. In 
future surveys, we will likely remove these latent priors 
to allow us to define upper limits on {Rs/Rp) as well as 
the usual (Ms/Mp). 

Although the photometry has been detrended and nor¬ 
malized using CoFiAM, we allow the out-of-transit base¬ 
line flux to be varied in our fits too. This term is lin¬ 
ear with respect to the transit model and thus can be 
minimized analytically at each Monte Carlo realization, 
which we implement here following the technique de¬ 
scribed in Kundurthy ct al. (2011). 

3.4. Vetting 

For each KOI studied, we obtain the Bayesian evidence 
and parameter posterior distributions for the planet-only 
model. Up, and also the planet-with-moon model, Rs- 
This information must be used to assess whether the KOI 
in question can be claimed to exhibit compelling evidence 
for the presence of one or more exomoons. In Kipping et 
al. (2013a), we first defined a list of standardized basic 
criteria that we use to make this assessment: 

B1 Improved evidence of the planet-with-moon fits at 
>4(7 confidence. 

B2 Planet-with-moon Bayesian evidences indicate 
both a mass and radius preference for the satellite. 


TABLE 2 

Priors adopted for the planet and moon parameters in our fits. 


Parameter 

Priori 

Planet Parameters 

p 

W{0,1} 

logio(p* [gcm-3]) 

W{-3,3} 

b 

W{0,2} 

P [days] 

W{P- 1,P-I- 1}* 

T [days] 

W{f- l,f + 1}* 

9? 

W{0,1} 

92" 

W{0,1} 

Moon Parameters 

(Rs/Rp) 


(Ms/Mp) 

W{0,1} 

asP 

W{2,7.897PA3}* 

COS i'g 

W{-1,3} 

Qs [rads] 

U{ — 7T, Tt} 

Ps [days] 

37(0.052, P/V3}* 

<I>S [rads] 

W{0, 27r} 


t U and ff denote a uniform and Jeffrey’s prior 
* P and r are the reported ephemeris parameters 


taken from the NASA Exoplanet Archive 
^ Defined in Kipping (2013) 

* Defined in Kipping et al. (2013b) 

B3 Parameter posteriors are physical. 

B4 Mass and radius of the moon converge away from 
zero. 

As discussed in §3.3, criterion B3 is enforced for all 
fits by virtue of our latent priors. This forces the mean 
density of the planet and the satellite to be within the 
extreme range discussed in Kipping ct al. (2012). It is 
therefore not explicitly tested. Criterion B1 remains un¬ 
changed in this work. Although Bayesian model selection 
is an important prerequisite for claiming a detection, it 
is insufficient in isolation since we adopt a Gaussian like¬ 
lihood function to compute the evidences, and real pho¬ 
tometry is never perfectly Gaussian. 

As of Kipping et al. (2013b), criteria B2 was modified. 
Formally, B2 requires running two extra photodynami- 
cal models in addition to 'Hsi which naturally slows down 
our survey speed considerably. In Kipping ct al. (2013b), 
we proposed focussing on the radius aspect of B2 alone, 
which we dub B2a, and testing the condition by allow¬ 
ing Rs to explore “negative radius” exomoons. Such 
signals correspond to inverted transits and we therefore 
require a positive radius moon as part of our vetting 
tests. Negative-radius moons have no physical meaning 
and are most likely indications of time-correlated noise 
structure. 

For criterion B4, the radius test is not useful for vetting 
when a latent prior exists on ps, such as done in this 
work, since zero-radius moons are forbidden. Instead, 
we focus on the mass component which remains unbiased 
and dub the criterion B4a. 

These modifications lead to the following basic detec¬ 
tion criteria: 

B1 Improved evidence of the planet-with-moon fits at 
> 4 (T confidence 

B2a Planet-with-moon evidences indicate a preference 
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for a positive radius moon (34.13% quantile of 
(Rs/Rp) is positive) 

B4a Mass of the moon converges away from zero (false 
alarm probability of a zero-mass moon is < 5% 
using the Lucy & Sweeney 1971 test). 

Additionally, Kipping ct al. (2013a) defined several 
“follow-up” criteria of which the most powerful is F2: 
“The predictive power of the moon model is superior (or 
at least equivalent) to that of a planet-only model”. This 
test first requires that we only fit ~ 75% of the transits 
in models Hp and Hs- We then a) extrapolate the max¬ 
imum a-posteriori realization of the two hypotheses and 
require Rs to yield a better likelihood and b) same as a), 
except we repeat for 10^ random draws from the poste¬ 
rior distributions of the two hypotheses and require Rs 
to be superior in 95% of the realizations. These respec¬ 
tively define criteria F2a and F2b. 

In practice, we only conduct F2b if F2a is first passed, 
to save resources. In some instances, the limited number 
of transits available meant we elected to fit 100% of the 
data, making F2 unavailable as a test. Otherwise, the 
^ 25% of ignored data is chosen to be near the median 
of the available LC and SC transits. 

F2 is a powerful tool in vetting exomoons but will be 
passed by KOIs exhibiting high signal-to-noise sinusoidal 
TTVs, since a moon-like TTV is also a sinusoid (Kipping 
2009a). KOI-0308.01 and KOI-269I.OI both show very 
strong evidence for TTVs (> 17 a) and nearly sinusoidal- 
like behavior, meaning they both passed F2a and F2b. 
As a more useful test in these two instances, we define the 
null model to be a finite-mass, zero-radius moon model 
Rs,p.o- This model is able to reproduce the sinusoidal 
TTVs but clearly a real exomoon case should be favored 
over a model negating the moon’s radius. 

3.5. Classification 

Naively, one might expect a KOI to either be classi¬ 
fied as a confirmed exomoon or a null result. In prac¬ 
tice, experience has revealed that spurious detections are 
common and are generally unsuitable for deriving upper 
limits. For this reason, it is necessary to classify null 
detections (“ND”) versus false-positives (“FP”). In our 
project, the upper limit we are interested in is the mass- 
ratio (Ms/Mp), since (Rs/Rp) is known to be both less 
sensitive and less reliable e.g. see Kipping ct al. (2013a,b, 
2015). Therefore, we define null detections as being cases 
where the posterior distribution of (Ms/Mp) does not 
converge away from zero. This is equivalent to criterion 
B4a, providing a quantitative method of classifying null 
detections. 

False-positives can be due to a variety of effects and 
in many cases we are unable to find an unambiguous 
explanation (except that we know an exomoon is 
unlikely to be responsible). These unknown false posi¬ 
tives are dubbed “FP-U”. In cases where instrumental 
effects are thought to be responsible, we dub the case 
“FP-I”. A recent example of this was discussed in detail 
for Kepler-90g.01 in Kipping ct al. (2015), where a 
pixel level effect was found to introduce an exomoon 
false-positive. Finally, if perturbations from another 
body are likely responsible, we dub the case “FP-P”. 
These are identified by assessing whether the TTVs 


show a significant signal, with > Aa confidence set as 
our threshold for classification as FP-P. 

4. RESULTS 

4.1. Exomoon Survey Results 
4.1.1. Overview 

Out of the 41 KOIs surveyed in this work, we find no 
compelling evidence for exomoons. Using our detection 
criteria discussed in §3.4, we find that only KOI-0072.01 
(Kepler-lOc) and KOI-3663.01 (PH-2b) pass these stan¬ 
dardized tests. However, in both cases we conclude that 
the moon candidates are false-positives, albeit for differ¬ 
ing reasons discussed shortly. 

In Table 3, we summarize the results of our survey 
for each KOI, providing the 68.3% credible intervals for 
(Ms/Mp) and (Rs/Rp) computed from our marginal¬ 
ized posteriors. We also detail which detection criteria 
were passed/failed and assign a classification (see §3.5) to 
each object. For some KOIs, two rows are shown due to 
two distinct modes in the parameter (Ms/Mp). In Fig¬ 
ure 2, we show the marginalized posterior of (Ms/Mp) 
for each KOI (again showing all modes if more than one 
present). 

4.1.2. Kepler-lOc 

For KOI-0072.02, also known as Kepler-lOc, we elected 
to use all of the Kepler data, rather than the usual 
75% strategy discussed in §3.4. This is because KOI- 
0072.02 is rare in providing us access to an additional 
vetting test. Specifically, a radial velocity solution ex¬ 
ists for this object (Dnmusque ct al. 2014), allowing us 
to compare the radial velocity derived planetary den¬ 
sity, pp, to that from our photodynamical fits. We find 
two modes for KOI-0072.02 for which the first is a clear 
null detection, failing criterion B4a. The second mode 
passes all of our basic detection criteria though, and 
here the pp test is of great value. The second mode 
requires pp = 17.1^^'3 gcm“^, which is highly incompat¬ 
ible with the radial velocity measurement of Dumusqiic 
et al. (2014) of pp = 7.1 ± l.Ogcm”^. This mode is 
therefore identified as a false-positive. Furthermore, as 
discussed in §4.2, KOI-0072.01 shows evidence for TTVs 
at the 5.0(7 level, meaning that the likely origin of this 
spurious detection is perturbations (i.e. FP-P). The con¬ 
sequences of TTVs in the Kepler-10 system are discussed 
in §4.2. 

4.1.3. PH-2b 

For KOI-3663.01, also known as PH-2b, all of the basic 
detection criteria are passed but no follow-up tests are 
possible due to the small number of transit epochs avail¬ 
able (just five). Plotting the maximum a-posteriori real¬ 
ization from model Rs reveals a solution driven by weak 
evidence for moon-like transits in the first four transit 
epochs (see Figure 3). As discussed later in §4.2, TTVs 
are not favored for this planet, since model Rp is pre¬ 
ferred over Rp^TTV at the 5 a level. As a result, the 
reality of the moon candidate rests solely on the reality 
of these moon-like transits. 

The third moon-like transit, present in QIO, occurs 
almost perfectly in sync with a data gap caused by 
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TABLE 3 

Summary of survey results. KOIs with multiple rows denote multiple modes. 


KOI 

(Ms/MpY 

(RsIRp) 

B1 

B2a 

B4a 

F2a 

F2b 

Classification 

Est. Ms [M®] 

KOI-0438.02 

< 0.96 

n 7 oc:+C.u:^« 
U.^OO_o 027 

/ 

X 

X 

N/A 

N/A 

ND 

< 5.07 

KOI-0518.03 

< 0.91 

+0.48t°-jl 

/ 

/ 

X 

N/A 

N/A 

ND 

< 5.51 

KOI-0854.01 

< 0.94 

+0.78i;;i® 

X 

/ 

X 

X 

- 

ND 

< 5.33 

KOI-0868.01 

(N O 

o o 
o q 
d d 
-h 1 

o 

cd 

-0.076«°}4 

X 

X 

/ 

N/A 

N/A 

FP-P 

- 

KOI-1361.01 

< 0.87 

+ 0 . 66 t°;i« 

X 

/ 

X 

X 

- 

ND 

< 5.36 

KOI-1431.01 

0 QQ 1 71+0 00024 

U.UUl i J-_o.00023 

— 0 141+0-015 

/ 

X 

/ 

N/A 

N/A 

FP-P 

- 

KOI-1783.02 

0 40+0-24 

u.zd^_q Q 37 

/ 

X 

/ 

N/A 

N/A 

FP-P 

- 

KOI-1830.02 

< 0.15 

4-0 445+0-085 

-hU.^^O_Q Qgg 

X 

/ 

X 

/ 

4 X 

ND 

< 1.14 

KOI-1876.01 

< 0.047 

—0 547+0-057 

U.O^/ _0 052 

X 

X 

X 

N/A 

N/A 

ND 

< 0.31 

KOI-2020.01 

< 0.0019 

-0.540tg°fg 

/ 

X 

X 

X 

- 

ND 

< 0.11 

KOI-2686.01 

n 1 ro+0.028 
U.iOC)_Q 

-0.286«°« 

/ 

X 

/ 

X 

- 

FP-P 

- 

- 

0 - 68 l°n 8 

U.ZDD_q Q 4 g 

/ 

X 

/ 

X 

- 

FP-P 

- 

KOI-2691.0lt 

0.63«Jg 

1 n 100-1-0.043 
-hU.l«4_Q Q 4 Q 

X 

/ 

/ 

X 

- 

FP-P 

- 

KOI-3263.01 

< 0.0020 

— 0 Q 73+0 017 

u.u/C)_o 018 

X 

X 

X 

N/A 

N/A 

ND 

< 0.041 

KOI-4005.01 

< 0.30 

-0.471«;«?4 

X 

X 

X 

/ 

X 

ND 

< 1.91 

KOI-4036.01 

< 0.048 

+0.581«;;f5 

X 

/ 

X 

X 

- 

ND 

< 0.26 

KOI-4054.01 

< 0.069 

n 71 o-l-0.085 
' ^'^-0.086 

X 

X 

X 

X 

- 

ND 

< 0.40 

KOI-0092.01 

U.ib 8 _Q gg 2 

4-0 Qn7+0 0t>t) 

-hU.CSU/_Q Ogg 

X 

/ 

/ 

X 

- 

FP-U 

- 

- 

< 0.017 

0 91 '7+0.037 
^•^1 ‘ -0.048 

X 

X 

X 

X 

- 

ND 

< 0.11 

KOI-0112.01 

< 0.26 

+0.306t°;M® 

X 

/ 

X 

/ 

X 

ND 

< 2.18 

KOI-0209.01 

< 0.010 

-0.30031°;°°^^ 

/ 

X 

X 

X 

- 

ND 

< 1.12 

KOI-0276.01 

< 0.65 

-0.26l°®t 

X 

X 

X 

X 

- 

ND 

< 4.31 

K01-0308.0lt 

u.»oz_Q 070 

0 oqo+ 0.020 
U.990_Q Q;Lg 

X 

X 

/ 

/ 

X 

FP-P 

- 

KOI-0374.01 

< 0.60 

-hU.ZUO_Q Q 3 g 

/ 

/ 

X 

X 

- 

ND 

< 5.03 

KOI-0398.01 

< 0.021 

0 1 01 -1-0.041 

-U. 1 » 1 _Q Q 3 g 

X 

X 

X 

X 

- 

ND 

< 1.47 

KOI-0458.01 

0 74+0.15 

b- '^- 0.20 

+0.857t°oT2l 

X 

/ 

/ 

X 

- 

FP-U 

- 

- 

< 0.88 

+ 0 . 66 l°Ji 

X 

/ 

X 

X 

- 

ND 

< 71.6 

KOI-0847.01 

< 0.024 

0 o'7o+0.023 

U.4/'5_q q25 

/ 

X 

X 

X 

- 

ND 

< 0.062 

KOI-1535.01 

< 0.012 

— 0 4 fiQ +0 049 
U.4Dy_Q Q 4 J 

X 

X 

X 

X 

- 

ND 

< 0.077 

KOI-1726.01 

< 0.17 

—0 41 

O.^l^_0.032 

X 

X 

X 

X 

- 

ND 

< 0.97 

KOI-1808.01 

0.35l°;ll 

4-0 705+0-110 

-hU.^DO_Q Qg 5 

/ 

/ 

/ 

X 

- 

FP-U 

- 

KOI-1871.01 

< 0.92 

-0.550t°;°®3 

X 

X 

X 

/ 

X 

ND 

< 5.36 

KOI-2065.01 

< 0.40 

+0.6691°;°^^ 

X 

/ 

X 

X 

- 

ND 

< 2.89 

KOI-2762.01 

< 0.78 

— 0 07‘^+0 055 

U.D/cJ_q Q 47 

X 

X 

X 

X 

- 

ND 

< 5.16 

KOI-5284.01 

< 0.84 

- 1-0 447+0-073 

-0.075 

X 

/ 

X 

N/A 

N/A 

ND 

< 2.51 

KOI-0072.02 

< 0.0025 

4-0 3i3+’-’-^Ji^ 

-|-U.O± 0 _Q 

/ 

/ 

X 

N/A 

N/A 

ND 

< 0.042LJ 

- 

n 1 90 + 0.043 

u.iZO_o 043 

4-0 423+0-036 
^'^•^^'^-0.031 

/ 

/ 

/ 

N/A 

N/A 

FP-P* 

- 

KOI-0245.01 

< 0.033 

4-0 329+0 016 
^^•'^^^-0.017 

/ 

/ 

X 

X 

- 

ND 

< 0.18 

KOI-0386.02 

< 0.0059 

4-0 434+0 027 

-|-U.^O^_Q Q 3 Q 

/ 

/ 

X 

X 

- 

ND 

< 0.046 

KOI-0518.02 

< 0.011 

-0.370tg;“f5 

X 

X 

X 

X 

- 

ND 

< 0.046 

- 

< 0.92 

+0.80lg;i3 

X 

/ 

X 

/ 

X 

ND 

< 4.02 

KOI-0722.01 

u.uo(_o 020 

4-0 584'*'®'^^'^ 

-hU.05^_Q Qgg 

X 

/ 

/ 

X 

- 

FP-U 

- 

KOI-1783.01 

< 0.030 

u.iuy_Q q29 

/ 

X 

X 

N/A 

N/A 

ND 

< 1.21 

KOI-2358.01 

< 0.015 

4-0 543+0-056 

-|-U.0^9_q Qg 4 

/ 

/ 

X 

X 

- 

ND 

< 0.090 

KOI-3663.01 

u.uu/a_o 0030 

4-0 121+0-013 

-hU.iZi_Q 

/ 

/ 

/ 

N/A 

N/A 

FP-I 

- 

KOI-3681.01 

< 0.0029 

4-0 121+0 006 
-hU.iZi_Q 231 

/ 

/ 

X 

N/A 

N/A 

ND 

< 1.84 


N/A: Test not available, since all observed transits were used in basic analysis. 

*: The 95.45% quantile is quoted for {M's/Mp) if B4a is not passed. 

^: Strong evidence for TTVs; null fit was 'Hs,ro rather than 'Hp 

Classified as a false-positive since pp from moon model is incompatible with independently measured pp 
Empirical planetary mass available from radial velocities 
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Fig. 2.— Histograms of the parameter {Mg/Mp) for the 41 KOIs surveyed in this work. Gray histograms are false-positive exomoon 
signals driven by some other perturber, colored ones are null detections. The estimated value of (Ms/Mp) is provided in the upper-right 
corner of each panel. 
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Fig. 3. — Transit light curves of all five epochs (each row) of PH-2b along with the maximum a-posteriori realization from model 'Hg 
(black line). Right panels show a zoom-in within ±0.8 d of the marked point of the left panel. Open diamond data points have zero weight 
in our fit as they have non-zero SAP QUALITY flags. The blue region denotes the 95.45% confidence interval of our model. The low 
significance of the model, plus coincidental location of two of the moon signals with instrumental effects, leads to us reject this candidate. 


points being ignored by CoFiAM, since the SAP QUAL¬ 
ITY flags are non-zero (recall that we do not include 
such points in our analysis, see §3.1). Specifically, from 
BKJDutc 926.93 to BKJDutc 927.34, there is a reac¬ 
tion wheel zero crossing. Without these points, there is 
in fact no data supporting the exomoon transit predicted 
by model 'Hs- Including the affected data show they are 
inconsistent with the model, as visible in Figure 3. 

For the second moon-like transit there is also a data 
gap obscuring the expected position of the moon’s ingress 
signal. Remarkably, as with the aforementioned case, 
these points were also ignored by CoFiAM due to another 
reaction wheel zero crossing event. Including the affected 
points shows that these data are consistent with model 
Tis- However, it should be cautioned that low amplitude 
non-periodic transit signals, such as this, are known to be 
reproducible by instrumental artifacts (e.g. see Kipping 
et al. 2015). 

A final blow to the exomoon hypothesis comes from dif¬ 
ferent detrendings of the same time series. For example. 


repeating CoFiAM on the PDCSAP data and comparing 
the maximum a-posteriori model from to the PDC¬ 
SAP detrended data reveals negligible improvement in 
the versus that of model 'Hp. Specifically, we find 
(Xp — Xs) = +1-7 for 3560 data points used in our anal¬ 
ysis. We repeated this test using a variant of the Trend 
Filtering Algorithm (TFA) designed for Kepler and de¬ 
scribed in Huang et al. (2013). Here, the results are even 
worse, with the planet model Up giving a superior ex¬ 
planation of the data with (xp — Xs) = —83.4. As a 
result, we consider this candidate to be an instrumental 
false-positive, reminiscent of the case of Kepler-90g.01 
discussed in Kipping et al. (2015). 

4.2. Transit Timing Variations 
4.2.1. Method 

We And that 12 of the 41 KOIs surveyed pass criterion 
B4a, suggestive of a perturbation acting on these ob¬ 
jects. To evaluate whether significant TTVs are present 
for these KOIs, we fit an additional planet model with 
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freely varied transit times, Hp^ttv- This model has the 
same 5 basic shape parameters for usual planets fits (see 
Table 2) with the global epoch term (r) and period [P) 
now fixed to the maximum a-posteriori result from model 
Tip. Additionally, for each of the N transits, an extra 
free parameter is included for that epoch’s time of tran¬ 
sit minimum, t^, leading to (5-1-A^) free parameters. For 
each KOI, we evaluate the significance of the fit using 
the usual Bayesian model selection method. The TTVs 
and the significance of the fit are shown in the right-hand 
panels of Figure 4. We also compute a periodogram us¬ 
ing the method described in Kipping ct al. (2013b) (es¬ 
sentially a Lomb-Scargle periodogram with a BlC-based 
definition of power) and plot the highest power signal 
over the TTV data. 

Since the speed of MultiNest scales poorly with di¬ 
mensionality (Handley et al. 2015), we find that KOIs 
with N > 17 were not practical to fit using model 
Tip, TTV- Instead, we fit each transit independently fixing 
the limb darkening coefficients to theoretical values inter¬ 
polated from the Claret & Blocmen (2011) tabulations 
and the period, P, as before, leaving 4 free parameters 
per fit. We refer to this model as Tip,i and show the 
results in Figure 5. Since the summed Bayesian evidence 
now accounts for excessive freedom, such as duration and 
depth changes, we instead define a signihcance for the 
TTV ht using the Bayesian Information Criterion be¬ 
tween the maximum periodogram peak signal and a null 
fit. These significances are not as statistically robust 
as the Bayesian evidence method, but serve as a useful 
proxy. 

4.2.2. Dynamically Active KOIs 

Seven of twelve systems inspected for TTVs show 
> 5a evidence for TTVs. Five of these show extremely 
strong evidence (> 17 a) for which there can be little 
doubt the KOIs are being perturbed; specifically KOI- 
0308.01, KOI-0868.01, KOl-1783.02, KOI-2686.01 and 
KOI-2691.01. Of these hve, only KOI-1783.02 is known 
to reside in a multiple transiting planet system and the 
TTV period is compatible with a perturbation from KOI- 
1783.01, although there is insufficient data to claim so 
unambiguously. The other four strongly perturbed KOIs 
are in single KOI systems, indicating at least four ad¬ 
ditional massive bodies not known to transit the parent 
star. 

Of the two marginal cases for TTVs, KOI-0072.02 and 
KOI-1431.01, the former is of particular interest being 
that this is Kepler-10c, the claimed solid Neptune-mass 
planet measured through radial velocities by Dnmusqne 
et al. (2014). Kepler-lOb is essentially decoupled from 
10c and is not expected to perturb 10c, suggesting an 
additional unknown massive body in the system. This 
object is likely to be of planetary mass to have avoided 
detection by radial velocities thus far. The period and 
mass of this object cannot be uniquely inferred from a 
single sinusoidal-like set of TTVs (Nesvorny & Vokrouh- 
licky 2014; Dock & Agol 2014). However, a near mean 
motion resonance configuration is most probable due to 
the fact that this configuration leads to the largest ampli¬ 
tude effects. Further, a planet on an exterior orbit would 
have a greater chance of having not been seen transiting 
to date, based on the system geometry. 

Whilst the evidence for TTVs in Kepler-lOc is not de¬ 


cisive, it is worth noting that the radial velocity peri¬ 
odogram (see Dumnsciuc ct al. 2014) shows numerous ad¬ 
ditional peaks to the frequencies of 10b and 10c. Notably, 
the next largest peak above the 0.1% false alarm prob¬ 
ability threshold occurs at 158.54 d, almost coincidental 
with a 7:2 ratio with 10c and potentially responsible for 
a near mean motion resonance interaction. We recom¬ 
mend additional radial velocities to identify the cause of 
the observed TTVs. 

4.3. Refined Planet Parameters 

An additional useful result from this work are param¬ 
eter posteriors derived for each of the planets surveyed. 
Since no exomoons are found, we use model Tip to esti¬ 
mate refined planetary parameters of each KOI studied. 
In the cases of the seven dynamically active KOIs, we use 
the TTV based models instead. The refined parameters 
are provided in Table 4. Figure 6 illustrates the differ¬ 
ences between the parameters derived in this work and 
those on the NASA Exoplanet Archive, which agree well 
in general. 

5. DISCUSSION 
5.1. Overview 

With the 41 KOIs surveyed in this work, the total num¬ 
ber of unique KOIs surveyed using Bayesian photody¬ 
namics is now 57. In every case, the analysis comes from 
the HEK project and were reported in Nesvorny ct al. 
(2012) and Kipping et al. (2013a,b, 2014a, 2015). In this 
survey, as with the previous smaller samples, we hnd no 
compelling evidence for an exomoon. 

Despite the sizable number of objects surveyed, we 
caution the community against attempting to derive 
constraints (the occurrence rate of exomoons) at this 
time, since some of our most interesting candidates re¬ 
main under detailed analysis and are yet to be presented. 
This is an inevitable consequence of the fact that pu¬ 
tative candidates require in-depth vetting, such as the 
case with Kepler-OOg.Ol where pixel-level analysis was 
required (Kipping ct al. 2015). 

5.2. Unknown False Positives 

Of the 41 KOIs studied, we find four instances of false- 
positives of unknown origin; KOI-0092.0I, KOT0458.01, 
KOI-0722.01 and KOI-I808.01. For the first three, cri¬ 
teria B2 and B4a are passed and for KOT1808.01 all of 
the basic criteria are satished. However, all four KOIs 
fail the follow-up test, F2a. This means that a Keplerian 
moon model provides a worse prediction (in a sense) 
of the ^ 25% of ignored transits than a vanilla planet- 
only model, indicating that a moon cannot explain the 
observations. 

When dealing with a handful of transits, quasi-periodic 
distortions to the transit profile, such as those due to 
spots (Bcky ct al. 2014), can be well fitted by the flexible 
exomoon model. However, since an exomoon is not the 
underlying cause, this model lacks any predictive power 
and thus should fail F2a. We therefore suggest that stel¬ 
lar activity is likely responsible for these four instances. 
As an example, we briefly discuss the most promising 
case of KOI-1808.01 which passes all of the basic criteria. 
Inspection of the posteriors reveals that the planet-moon 
separation is required to be tight at planetary 
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Fig. 4.— Right column shows the TTVs computed for 6 of the 12 KOIs which pass criteria B4a. These 6 KOIs have TV < 17 epochs 
enabling us to conduct the global fit, 'Hpttv- The left panel shows the corresponding Lomb-Scargle periodogram with a BIC-based 
definition of power. The favored trial is plotted over the TTV data on the right panel, with the significance denoted in the corner. Negative 
significances indicate that the TTV model is disfavored. 


radii leading to no clean, separated moon transits in the 
best fit. Instead, the signal is driven by distortions to 
the planetary transit event, consistent with the effects of 
spots. 

5.3. TTV Systems 

Although it is not the objective of our work, we find 
seven systems classed as FP-P, that is exomoon false- 
positives due to perturbations from another body. These 
cases display strong evidence for TTVs (> 5 a) and ac¬ 
count for 17% of the KOIs surveyed. Curiously, five of 
these instances occur for the 16 targets selected through 
TSA, which has a coincidental probability of occurring 


of just 7.7%. Since the TSA targets were selected to be 
in the habitable zone, they are more biased towards long 
period KOIs, perhaps suggesting that TTVs are more 
common for such worlds. Whatever the cause, based on 
the 41 KOIs studied here, habitable zone TSA planetary 
candidates have an a-priori chance of displaying TTVs 
of (31 ± 12) %. 

This rate is substantially higher than that reported by 
Mazch ct al. (2013), who find 130 KOIs with significant 
TTVs out of a total sample of 1897 i.e. (6.9 ± 0.6) %. 
This is incompatible with the rate seen for the TSAs be¬ 
yond the inner edge of the habitable zone (as defined 
by Zsoni et al. 2013) at a level of 3.2 tr, meaning co- 
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Fig. 5. — Right column shows the TTVs computed for 6 of the 12 KOIs which pass criteria B4a. These 6 KOIs have Ai > 17 epochs 
and thus have too many free parameters for the global fit, Hp ttv-. leading us to use a model of N independent fits, 'Hp i- The left panel 
shows the corresponding Lomb-Scargle periodogram with a BIC-based definition of power. The favored trial is plotted over the TTV data 
on the right panel, with the significance denoted in the corner. 


incidental sample selection is highly improbable to be 
responsible. However, the sample studied by Mazeh et 
al. (2013) casts a much wider net than the carefully se¬ 
lected targets of TSA. Specifically, this rate is for the 
2321 KOIs in the Borucki et al. (2011, 2012) catalogs 
with a SNR> 2.5 (giving 1960 KOIs) exhibiting 7 or 
more transits in the first 12 quarters of data (giving 
1897 KOIs). The latter filter excludes long period plan¬ 
ets (P < 150 days) and moreover the sample is heavily 
biased towards short period events due the higher detec¬ 
tion probability (^ p-5/3. Gaudi 2006). 

The tension between these two rates could be resolved 
by a period dependency for the occurrence rate of de¬ 


tectable TTVs. Our TSA sample is designed to pick 
gaseous planets in the habitable zone or further out, im¬ 
plying that habitable zone gaseous planets are frequently 
accompanied by additional planets, likely in period com- 
mensurabilities. These hidden planets probably reside 
in exterior orbits, in order to explain having not been 
seen to transit (since multi-planet systems are typically 
coplanar; Fabryckv’ et al. 2014). However, we caution 
that our TSA method uses a complex series of criteria, 
such as invoking tidal evolution models, which may make 
them not representative of habitable zone gaseous plan¬ 
ets. The occurrence rate of dynamically active habitable 
zone or cooler gaseous planets is therefore an interesting 
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TABLE 4 

Planet properties as derived using model Hp ■ 


KOI 

P 

loglo(P*[kgm 3]) 

b 

P[d] 

T [BJKD] 

91 

12 

KOI-0438.02 

U.UCii^/_0.00089 

4.158i;j;f,« 

0 96+'i'^^ 
U.OD_o.25 

52.6615991?-°????; 
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0 25+0.36 
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KOI-0518.03 

n 0341 c;+0 00229 

U.U04HD_o. 00088 

3-44l°i° 
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0.48l?-°; 

0 29+0-31 
6-^'^-0.15 

KOI-0854.01 

n rioQQ+ 0.0020 
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3 fA9+0-14 
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0 46"*”6-61 
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56.056081?;????? 

873.83701?-???? 

0 97+0.37 
'-'•^^-0.19 

0 Q-i +0.38 
6-'51-0.22 

KOI-0868.0lt 

0 151 fi+00024 

U.i01D_o.ooi4 

9 o'rrj+O.OlO 

0.4/U_o.oi4 
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680.404401?;????° 

0 oc::+0.30 

n 9 Q + 0.36 
6-^0-0.20 

KOI-1361.01 

U.U34Q0_Q OQQgg 

o.uzu_o_284 

n q7+0.30 
^•'^^-0.25 

“iQ 07709+0.00012 
oy.o/ '^^-0.00013 

869.72081?-??;; 

0 C9+0.29 

y-O^-0.23 
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6 -^o_o.i7 
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0 0774 + 0-0026 

U.U^ <^_0.0014 
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*^•"^^^-0.032 
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34^1 1 5Qr)2'*”6‘66035 
CS40.10yU4_Q QQQ 35 
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^• 6 ^- 0.12 

0 20+0-38 

6 -^y- 0.21 

KOI-1783.02t 
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q 1 Q+0.52 
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0 oq+0.36 
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0 55+0.29 
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0.401?;?; 
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b.bU_o.28 

KOI-2020.01 

0 QQ7Q+0-0039 
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'^•^^-0.40 
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0 4Q+0-31 
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0 6tS+6-^4 
b.bO_o.3o 

KOI-3263.01 

0 zlS + 0-^'^ 

0-4o_o.29 

4 924+0-064 
4.444_o 055 

1 90+0.39 

i. 20 _o 37 

'7f\ Q7Qqqq+0.000046 
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KOI-4005.01 
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KOI-4054.01 
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878.34721?-°??? 

0 1 g+0.29 
^•-*-^-0.13 

n OO+0.40 
6-Oo o.24 

KOI-0092.01 

u.uz^y_0 0018 

9 Qn+^-^6 

'9-oU_0 30 

0 71+^-i^ 

^1-0.42 

65.7045651?;?????? 

70X4^000+11^110032 

t y^.^OOU9_Q QQQ32 

0 C2+6-^^ 

'-'•6^-0.16 

0 17+0-^J- 

6 -l'- 0.11 

KOI-0112.01 
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U.U 2 » 2 _o 0015 

^•'^^-0.13 

0 784+6-645 
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- 

- 

KOI-0374.01 
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'^-0.27 

0 61''"6-16 
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755.03681?-°?;? 

0 45+^-^^ 
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KOI-0398.01 

u.uy^o_o 0018 

o.^oi_o 045 
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Cl 8/lfi8^8+6-000028 
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0 75+0.17 

' 6 - 0.21 
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U.OOU_Q Qgg 

KOI-0458.01 
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^•''^- 0.12 
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00 ./ 10 UZ/_Q QQQQgg 

906.414421?;????? 
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KOI-0847.01 

0.0568t°““3 

2 664"'"^'^^^ 

Z.OO^_OJ^55 

0 41+0.19 
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80.871991?;???;? 

850.872361?;????? 

0 41+6-^i 
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6-^^-0.17 

KOI-1535.01 
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'=>•^^-0.58 
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- 

- 
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O.Uoo_o 012 

0 5950+6 OO84 
U.JOOU_Q QQgg 
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^.yo_o.45 
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'HpTTV used instead due to significant TTVs present 
*: Plpj used instead due to significant TTVs present and >17 transit epochs 


question to pursue in future work. 

5.4. Empirical Sensitivity to Exomoons 

Several previous works have attempted to estimate the 
expected sensitivity of Kepler and other telescopes to ex¬ 
omoons (e.g. Kipping et al. 2009; Simon et al. 2009). 
None of these works have gauged the sensitivity of full 
photodynamics, rather focusing on just one or two of the 
multitude of effects by which an exomoon may betray 
its presence. Further, real noise structure and complex 
inter-parameter degeneracies have not been realistically 
modeled. The most robust estimate of our photomet¬ 
ric sensitivity to exomoons comes from examining the 
empirical upper limits derived from a reasonably large 
sample, such as those from the HEK project surveys. 


Of the 57 unique KOIs studied, we find 46 null detec¬ 
tions which are suitable for deriving upper limits from. 
This indicates that we obtain null detections for (81±5)% 
of cases studied, with the others being false-positives. 
From the null detections, we find a diverse range of sen¬ 
sitivities in {Ms/Mp), ranging from 1.9 x 10“^ (95.45% 
upper limit of KOI-2020.01) to nearly unity, as shown 
in Figure 7 (for reference, the Earth-Moon mass-ratio is 
12.3 X 10“^). This diversity is likely due to a complex 
combination of factors, such as the target magnitude, 
planetary transit depth and duration, correlated noise 
structure. Hill sphere size, latent TTVs from perturb¬ 
ing planets, photometric cadence and number of transits. 
Regardless, we are able to infer that, assuming our sam¬ 
ple is typical, the null detections are sensitive to Earth- 
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(a/H*) 


Fig. 6. — Comparison of the NASA Exoplanet Archive listed 
transit parameters (triangles) versus the values derived in this work 
(circles). Specific KOIs are connected by a line and have the same 
color coding. The plotted points represent the medians of the pa¬ 
rameter posteriors. 


Moon mass-ratios or better for ~ 15% of cases and to 
Pluto-Charon mass-ratios for ~ 47% (based on 95.45% 
confidences). 

The relative mass limits of (Ms/Mp) may be con¬ 
verted to absolute masses limits in Ms by estimating 
Mp. For cases such as Kepler-10c, the physical mass is 
known (Dumusquc et al. 2014), and here our upper limit 
corresponds to 1.7 Ganymede masses (at 95.45% confi¬ 
dence), demonstrating our ability to find moons akin to 
those in the Solar System. In the other cases, we must 
assign a mass based on the observed radius, Rp. For 
self-consistency, we use the same Rp-Mp relation used 
for TSA, namely that discussed in Kipping et al. (2014a), 
which is predominantly an empirical relationship. As 
shown in the right panel of Figure 7, this allows to esti¬ 
mate that we are sensitive to Earth-like moons or smaller 
for ~ 40% of the null cases. The lowest mass exomoon 
that can be considered capable of supporting an Earth¬ 
like atmosphere is thought to be at 0.3 (Raymond et 
al. 2007), for which we find sensitivity for ~ 30% of null 
cases. 

In summary then, based on empirical sensitivity limits, 
we show for the first time that the HEK project is sensi¬ 
tive to even the smallest moons capable of being Earth¬ 
like for I in 4 cases (after accounting for false-positives). 
In terms of planet-mass ratios, we find even that the 
Earth-Moon mass-ratio is detectable for 1 in 8 of cases, 
posing a challenge but not an insurmountable barrier. 
Mass ratios of ^ 10“"^, such as that of the Galilean satel¬ 
lites, have never been achieved. However, if Galilean-like 
satellites reside around lower-mass planets than Jupiter, 
of order ~ 20 M 0 , then we do find sensitivity, as demon¬ 
strated by the limit of 1.7 Ganymede masses achieved for 
Kepler-1 Oc. 


5.5. Computational Demands 

We briefly discuss the computational demands of the 
HEK survey to date. Previous surveys generally used 
less than the full Kepler time series (Q0-Q17) due to the 
timeline of the data releases. In contrast, all 41 KOIs sur¬ 
veyed here used all available quarters, providing a fairer 
assessment of computational requirements for future sur¬ 
veys. We recorded the computational time used for each 




Fig. 7. — Upper: Smoothed empirical sensitivity of the HEK 
survey to the satellite-to-planet mass-ratio. Based on 46 null re¬ 
sults (thin gray line) drawn from the 57 unique KOIs studied. HEK 
is sensitive to Earth-Moon mass-ratios for ~ 15% of null cases and 
to Pluto-Charon mass-ratios for ~ 47% (based on 95.45% confi¬ 
dences). Lower: Same as top, except sensitivity quoted in abso¬ 
lute masses via estimated planetary masses. HEK is sensitive to 
0.3 M0 mini-Earths for ~ 30% of cases and to Earths for ~ 40% 
(based on 95.45% confidences). 

model fit Hs, in order to better understand our require¬ 
ments. A histogram of the CPU times is shown in Fig¬ 
ure 8, with a mean time of 33,000 hours per KOI. Note 
that these times do not include the planet-only fits (neg¬ 
ligible times) and TTV fits (sizable but less demanding). 
For all fits, the fits were conducted on AMD Opteron 
62xx/63xx CPUs. The spread in times appears mostly 
linked to the number of transit epochs fitted, with larger 
time series naturally requiring more overhead. With cur¬ 
rent in-house server capacity, we estimate that we are 
able to survey at least 100 KOIs per year. 

5.6. The Effects of Correlated Noise 

Despite the groundbreaking precision achieved by Ke¬ 
pler (Christiansen et al. 2012), when one seeks low 
signal-to-noise events then the effects of time-correlated 
noise become important. This is particularly salient 
for exomoons, since the transits do not follow a simple 
ephemeris which can be used to check for repeatability 
(Kipping 2011). This allows the model great flexibility 
in explaining any slight random changes in flux, leading 
to spurious detections (e.g. Kipping ct al. 2015). 

One strategy to guard against correlated noise is to 
use a likelihood function more suited for such structure. 
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Fig. 8.— Histogram of the number of CPU hours used for model 
1-is oi the 41 KOIs surveyed in this work. The distribution has a 
mode of 24, 000 CPU hours and a mean of 33, 000 CPU hours. 


for example using wavelets or Gaussian processes (ped¬ 
agogical examples of these methods applied to photom¬ 
etry can be found in Carter & Winn 2009 and Gibson 
et al. 2012, respectively). However, these methods have 
so much flexibility than they may fit-out the very low 
amplitude signals we seek. An alternative strategy, as 
employed in this work, is to use vetting criteria, such as 
the double likelihood test F2 and insights from our photo- 
dynamical fits. Whilst maintaining the use of a Gaussian 
likelihood leads to more false-positives, it at least allows 
us see the full array of candidate signals and to assess 
them accordingly. Since the use of Gaussian likelihood 
functions is so prevalent within the community, we here 
discuss the rate of false-positives induced by correlated 
noise structure. 

Correlated noise directly affects the time series pho¬ 
tometry, potentially mimicking transit-like signals and 
leading to spurious detections of {Rs/Rp)- Whilst cor¬ 
related noise can also affect the TTVs by distorting the 
planetary transit profile (Mazeh ct al. 2015), such transit 
distortions are usually fitted out by the radius effect of 
our moon model, as seen in the example of Figure 9. In 
what follows, we focus on the effect of correlated noise 
on moon transit signals, i.e. parameter {Rs/Rp). In our 
survey, a non-zero {Rs/Rp) moon fit can be caused by 
two distinct reasons: 

A] The latent prior of ps < 27.9 gcm“^ (see §3) tech¬ 
nically forbids {Rs/Rp) = 0, thereby forcing gen¬ 
uinely null {Rs/Rp) cases to be either positive or 
negative. 

B] Time correlated noise may drive the fits towards a 
false-positive likelihood minimum, since we adopt 
a Gaussian likelihood function. 

It is of particular interest to understand how frequently 
case B occurs, which severely affects all methods seek¬ 
ing evidence of exomoons via their transits. In hind¬ 
sight, not imposing the latent ps prior would be the 
most effective way to measure the occurrence of spurious 
{Rs/Rp) detections due to correlated noise. However, 
we discuss here an alternative approach using the infor¬ 
mation in hand. The mean satellite density is computed 
using Equation 5 of Kipping (2010b): 





BKJDutc 1554.9 


Fig. 9.— Example of a “negative radius” (inverted transits) 
exomoon solution, which occurs for around half of all KOIs sur¬ 
veyed. The transit light curves used for fitting KOI-4054.01 are 
shown, along with the maximum a-posteriori realization of our 
planet+moon model (black line) and the 95.45% credible inter¬ 
vals (blue shaded region). This model is favored over a planet-only 
at 3.5(J with {Rg/Rp) = —0.71di0.09, driven by correlated noise. 


PS = Pp 


{Ms/Mp) 

\Rs/Rp?' 


( 1 ) 


Since we impose the condition that ps < 27.9gcm“^, 
then plotting the joint posterior of {pp{Ms/Mp)) and 
\Rs/Rp\^ will have an excluded region defined by a 
straight line passing through the origin of gradient 
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27.9gcm“^. We inspected each of the joint posteriors 
as described, identifying cases where the posterior was 
bunched up against this boundary, indicating the latent 
prior was likely driving a non-zero (Rs/Rp) solution (i.e. 
case A). Of the 41 KOIs surveyed, we find that 22 show 
no evidence for being affected by the latent ps prior and 
therefore the non-zero (Rs/Rp) solution is reasoned, by 
deduction, to be caused by time correlated noise. 

This estimate provides a handle as to the prevalence 
of correlated noise-induced spurious moon-like transits, 
indicating that (54 ± 8)% of KOIs surveyed are affected. 
In future work, we intend to relax the latent priors and 
conduct a more in-depth analysis of this question. Of 
the 22 identified cases, 8 yield negative radius moon solu¬ 
tions, compatible with the expectation of 11±2.3, assum¬ 
ing that correlated noise induces negative radius transits 
with a p-value of 0.5. 

If we had sought evidence for exomoons by seeking 
exomoon radius signals alone, and had not used a time- 
correlated noise likelihood, we would have erroneously 
inferred evidence for exomoons in around half of these 
cases, which is over one quarter (~ 26%) of all KOIs sur¬ 
veyed. Thanks to photodynamics and our vetting criteria 
we were able to identify these spurious cases. 

Although our survey strategy is robust against these 
time-correlated noise-induced moon-like transits, other 
exomoon search methods that do not use photodynam¬ 
ics may struggle to identify these frequent occurrences. 


We highlight these cases as a service to the community, 
to emphasize that exomoons live in the regime where cor¬ 
related noise is present and one must employ methods to 
guard against it when seeking such signals. 
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